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Abstract 

An approach which closely maintains the non-dissipative nature of classical fourth or higher- 
order spatial differencing away from shock waves and steep gradient regions while being capable 
of accurately capturing discontinuities, steep gradient and fine scale turbulent structures in a stable 
and efficient manner is described. The approach is a generalization of the method of Gustafsson 
and Olsson and the artificial compression method (ACM) of Harten. Spatially non-dissipative 
fourth or higher-order compact and non-compact spatial differencings are used as the base schemes. 
Instead of applying a scalar filter as in Gustafsson and Olsson, an ACM like term is used to 
signal the appropriate amount of second or third-order TVD or ENO types of characteristic based 
numerical dissipation. This term acts as a characteristic filter to minimize numerical dissipation for 
the overall scheme. For time-accurate computations, time discretizations with low dissipation are 
used. Numerical experiments on 2-D vortical flows, vortex-shock interactions and compressible 
spatially and temporally evolving mixing layers showed that the proposed schemes have the 
desired property with only a 10% increase in operations count over standard second-order TVD 
schemes. Aside from the ability to accurately capture shock-turbulence interaction flows, this 
approach is also capable of accurately preserving vortex convection. Higher accuracy is achieved 
with fewer grid points when compared to that of standard second-order TVD or ENO schemes. 
To demonstrate the applicability of these schemes in sustaining turbulence where shock waves are 
absent, a simulation of 3-D compressible turbulent channel flow in a small domain is conducted. 
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I. Introduction 

Modern shock-capturing methods such as total variation diminishing (TVD) or essentially 
non-oscillatory (ENO) types of schemes that are higher than third-order accurate are usually 
CPU intensive, involve large grid stencils and require special treatment near boundary points. In 
spite of their high-resolution shock-capturing capability, these schemes often exhibit undesirable 
amplitude and/or phase errors for vortical and turbulent convection flows and complex wave 
propagation phenomena. See NASA Conference Publication 3300, May, 1995, and Sandham & 
Yee (1998) and references cited therein for some discussion. To compensate for the somewhat 
ad hoc ways of utilizing TVD or ENO schemes for compressible viscous flows, Toro (1991) 
proposed a viscous flux limiter approach to deal with scalar mixed hyperbolic-parabolic problems. 
Systemic extension of Toro’s idea to a system of equations containing other than a single scalar 
viscosity term remains a challenging area of research. The objective of this paper is to propose a 
compromise between the above two approaches while maintaining an efficient way of obtaining 
fourth or higher-order accuracy almost everywhere without using higher than third-order TVD or 
ENO schemes. Hereafter we refer to ‘ ‘high order schemes’ ’ as schemes with spatial accuracy that 
is greater than three away from shocks and steep gradient regions. 

Accurate and efficient direct numerical simulation (DNS) of turbulence in the presence of 
shock waves represents a significant challenge for numerical methods. A numerical scheme for 
DNS of shock-turbulence interactions of high speed compressible flows would ideally not be 
significantly more expensive than the standard fourth or sixth-order compact or non-compact 
central differencing scheme. It should be possible to resolve all scales down to the order of the 
Kolmogorov scales of turbulence accurately and efficiently, while at the same time being able to 
capture steep gradients occurring at much smaller scales. Appropriate numerical schemes should 
not interfere with the turbulence mechanisms resulting directly from the governing equations. See 
Sandham and Yee (1998) and references cited therein for a discussion. 

Gustafsson and Olsson (1995) developed stable high order centered schemes with stable 
numerical boundary condition treatments. For problems containing shocks, they used a scalar 
shock-capturing filter. Such schemes have advantages over higher-order ENO schemes which 
require very large grid stencils even for modest orders of accuracy. (For example, a seven-point 
grid stencil is required for a second-order ENO scheme.) In this paper we propose to use the 
narrow grid stencil of high order classical spatial differencing as base schemes. Low order TVD 
or ENO dissipation in conjunction with the Harten artificial compression method (ACM) switch 
(Harten, 1978) are used as characteristic filters. The ACM procedure is similar to Harten (1978) 
but applied in a slightly different context. The final grid stencil of these schemes, for example, 
is five if second-order TVD schemes are used as filters and seven if second-order ENO schemes 
are used as filters for a fourth-order base scheme. Numerical boundary condition treatment is 
simple and can be the same as for the existing base and filter schemes. Here, we propose to use 
filter operators that have similar grid stencil widths as the base scheme for efficiency and ease of 
numerical boundary treatment. Higher than third-order filter operators are of course applicable, 
but they are more CPU intensive and require special treatment near boundary points for stability 
and accuracy. On one hand, it would defeat the purpose of achieving efficiency. On the other 
hand, near shocks and shears, the resolution of higher than third-order TVD or ENO schemes 
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is comparable to their lower-order cousin except with a slight gain in resolution in regions near 
steep gradients and smooth flows. If lower order filters are able to minimize but at the same time 
provide the proper amount of numerical dissipation away from shocks and shears to stabilize the 
non-dissipative nature of the high order base scheme, we would achieve similar resolution with 
improved efficiency. This is the philosophy used to design the schemes. Our approach is aimed 
mainly at problems containing vortex convections, shock, shear, and turbulence interactions. As 
illustrated in later sections, these types of characteristic TVD (and ENO) filters can even improve 
fine scale flow structure when applied to existing methods of Harten (1983, 1984), Yee (1985, 
1987) and Yee (1997). 

Section II describes the numerical algorithm. Section III illustrates the performance of these 
algorithms for a variety of unsteady flows where most conventional methods exhibit difficulty 
in obtaining low dissipative solutions in an efficient and stable manner. The first problem is a 
stationary vortex evolution. The second problem is a convecting vortex. In the third problem, 
a vortex pairing in a time-developing mixing layer, shock waves form around the vortices. In 
the fourth problem, a shock wave impinging on a spatially evolving mixing layer, the evolving 
vortices must pass through a shock wave, which in turn is deformed by the vortex passage. To 
demonstrate the applicability of these schemes in sustaining turbulence where shock waves are 
absent, a simulation of compressible turbulent channel flow in a small domain is carried out. For 
problems 3 and 4, the detailed physics and extensive evaluation of the proposed scheme were 
reported in a separate paper by Sandham and Yee (1998). Here, only certain aspects of the 
performance of these schemes for the two problems are described. The study of the performance 
of this approach for time-marching to the steady-state numerical solutions is in progress. 


II. High Order Shock-Capturing Schemes Using Characteristic Filters 

For simplicity of presentation, the discussion will concentrate on the convection part of the 
Navier-Stokes equations. Analogous order of accuracy of spatial discretizations for the viscous 
terms will be briefly described at the end of this section. 

In vector notation the 2-D compressible time-dependent Euler Equations in conservation form 
for an equilibrium real gas can be written as 

Ut + F m + Gy = 0, (2.1a) 

where U t = ^ , F x = §~~ and G y = and the U, F, G, vectors given by 


u = 

i 

; F = 

pu 

pu 2 + p 

; G = 

pv 

puv 

2 i 


j 

<u 


puv 

. eu + pu . 


P V "1"P 
. ev + pv . 


T 

The dependent variable U is the vector of conservative variables, and (p, it, i>,p) is the vector of 
primitive variables. Here p is the density, u and v are the velocity components, pu and pv are the 
x- and ^-components of the momentum per unit volume, p is the pressure, e = p[e + ( u 2 + v 2 )/2] 
is the total energy per unit volume, and e is the specific internal energy. 
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For a thermally perfect gas, the equation of state is 

P = PRT, (2.2) 

where R is the specific gas constant, and T is the temperature with e = e(T). For constant specific 
heats (calorically perfect gas) 


e = c v T , (2.3) 

where c v is the specific heat at constant volume. 

The eigenvalues associated with the flux Jacobian matrices of F and G are (u,u,u ± c) and 
(v,v,v ± c), where c is the sound speed. The two u,u and v,v characteristics are linearly 
degenerate. Hereafter, we refer to the fields associated with the u ± c and v ± c characteristics 
as the nonlinear fields and the fields associated with the and v f v characteristics as the linear 
fields. 

The basic idea of these shock-capturing schemes consists of two steps. The first step is the 
high order spatial and temporal base scheme. Many of standard high order non-dissipative or low 
dissipative base schemes fits in the present framework. The second step is the appropriate filter 
for stability, shocks, contact discontinuities, and fine scale flow structure capturing. Many of the 
TVD and ENO dissipations, after a minor modification, are suitable candidates as filters. 


2.1. The Base Schemes 

In this paper, only the method of lines approach is considered. We divide the discussion of the 
base schemes into time and spatial discretizations. The filter step either does not directly involve 
the time discretizations or it uses the same time discretizations as the base scheme, depending on 
the types of temporal schemes. 


2.1.1. Time Discretizations 

Third or higher-order linear multistep methods (LMMs) (Gear, 1971; Lambert 1973) usually 
involve more than three time levels and initial starting schemes are required. For stiff problems, 
stiffly stable implicit methods are desirable, especially for time-marching to steady-state numerical 
solutions. Examples of explicit LMMs are explicit Euler and Adams-Bashforth. Examples of 
implicit LMMs are backward Euler, trapezoidal rule, and three-point backward differentiation. 
For non-stiff or moderately stiff multidimensional problems, one of the easiest procedures for 
obtaining higher than second-order time discretizations is the Runge-Kutta method. There are 
many variants of the Runge-Kutta method in the literature. See Lambert (1973), Butcher (1987), 
Carpenter & Kennedy (1994) and Gottlieb & Shu (1996) for details. Let 


U t = L(U) jjh (2.4) 

be the semi-discrete form of (2.1), where L is the spatial discretization operator for (— F x - G y ) 
to be discussed in the next section. If viscous terms are present, L includes the viscous spatial 
discretizations. Here Uj,k is a discrete approximation of U at x = j Ax and y = kAy, where Ax 
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and Ay are the grid spacing in the *- and y-directions and j and k are the corresponding spatial 
indices. For simplicity of discussion, uniform Cartesian grids are assumed. Generalization of the 
method to nonuniform grids with analytical coordinate transformation is straightforward. 

The fourth-order classical Runge-Kutta method takes the form 


Jfej = L(U n ) 

k 2 = I,(v n + 

= L(u n + 
k i = L(u n + Atk jj 

— U n + — f&i + 2 h 2 + 2k 2 -+■ ki . 


(2.5) 


Shu’s third-order Runge-Kutta (Shu, 1989) form that is compatible with TVD, TVB (total variation 
bounded) and ENO schemes takes the form 


U w = IT" + A tL(U n ) 

XjW = Ip- + ipd) + 1 a tL(U w ) 

4 4 4 

U n+1 = -U n + - U w + tL{U w ). (2.6) 

3 3 3 

The order of the temporal discretization might not be the key measure of the choice of temporal 
method. At times, one may be mainly interested in the phase error of the solution. Schemes which 
have higher-order accurate phase error might have lower order when measured in the standard L 2 
norm. For hyperbolic and wave-like problems, one usually desires the accuracy in time and space 
to be equal. Another consideration is that the combined spatial and temporal discretizations might 
pose a very stringent time step constraint for the overall scheme. In addition, the proper choice of 
time discretization that is compatible with a chosen spatial discretization is crucial in achieving 
low phase and amplitude errors for time-accurate computations. This subject is ongoing research. 
For all the model test problems considered in this paper, the classical fourth-order Runge-Kutta 
method appears to work well. 

2.1.2. Spatial Discretizations for the Convection Terms 

Denoting F^h as the discrete approximation of the convection flux F at (j Ax y kAy), samples 
of the high order base scheme for F x (similarly for G y ) can be of the following four types. 

Central Differencing s: (fourth and sixth-order) 


0 


F m 


12A* 


( F > 


i+2,fc — &Fj+i t k + — Fj- 2 




(2.7) 


F» 


60Az 


^-Pj+ 3 ,fc — 9F i+2 ,* + 45Fj + i,* — 45Fj*_ lt fc -f 9Fj_ 2 ,j, — Fj- 3 ,*^ . 


( 2 . 8 ) 


Compact Central Differencings: (fourth and sixth-order, Hirsh (1975), Ciment & Leventhal 


(1975) and Lele(1992)) 

F **~hc{ A ‘ lBmF ). k (29a) 

where for a fourth-order approximation 

(A m F)j,h = g (Fj + i t h + 4Fj,fc + F,*_ (2.9b) 
(*.*)*» = J (-Fi+x,* - Pi- i.i.) . (2.9c) 

and for a sixth-order approximation 

(A,F) j<k = | (fi+j,* + 3fj,» + f,-,.*) , (2.9d) 

(B.F) jtk = i + 28F,+i,j, - 28 F,- X jk ~ ^- 2 ,*) • (2-9e) 

Predictor- Corrector Differencings : (Fourth and Sixth-order) 


Predictor: ^ - 7F Jf fc + 8Fj-i,* - F/- 2 ,fc^ , (2.10a) 

Corrector: -^(7^ - 8F J+ i,/, + F i+2 , fc ), (2.10b) 

Predictor: ^ - 37 Fj,* + 45Fj~i,* - 9 Fj~ 2 ,h + Fj- 3fh ^ , (2.11a) 

Corrector: ^37Fj, fc - 45Fj +lik + 9F J+2 , fc - F J+ j |fc ^ (2.11b) 
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New forms of the upwind biased predictor-corrector methods including compact formulations 
recently developed by Hixon and Turkel (1998) are also applicable. Interested readers should 
refer to their paper for the various upwind biased predictor-corrector formulae. The choice of the 
time integrators for these types of predictor-corrector methods is more limited. For example, if 
second-order time accuracy is desired, then (2.10) and (2.1 1) in conjunction with the appropriate 
second-order Runge-Kutta method are analogous to the familiar 2-4 and 2-6 MacCormack schemes 
developed by Gottlieb & Turkel (1976) and Bayliss et al. (1985). Here the first number refers 
to the order of accuracy for the time discretization and the second number refers to the order 
of accuracy for the spatial discretization. However, in this case one achieves the second-order 
time accuracy without dimensional splitting of the Strang type (Strang, 1968). For higher than 
second-order time discretizations, only certain even stage Runge-Kutta methods are applicable. 
For compatible fourth-order Runge-Kutta time discretizations, see Hixon and Turkel ( 1 998) for 
possible formulae. For example, the classical fourth-order Runge-Kutta is applicable provided one 
applies the predictor and the corrector step twice for the four stages; i.e., the predictor step for the 
first and third stages and the corrector step for the second and fourth stages. 

The SHOEC Differencing s: (Gerritsen, 1996) The split high order entropy conserving scheme 
(SHOEC) of Gerritsen (1996) extends the summation by parts and entropy splitting idea of Olsson 
(1995) to the 2-D Euler equations for an ideal gas. It is based on the entropy splitting of the 
convection flux using Harten’s symmetrized form via entropy variables (Harten, 1981). Using the 
entropy variable transformation W = W(U), one splits, for example, F m 

F ■ = + Tb FwW ‘ (2 - 12a) 

with (3^-1 and Fw = The vector W is chosen such that both F(U(W)) and U(W) are 
homogeneous functions of the appropriate order /3. For the perfect gas 2-D Euler equations W, 
and Fw and Gw are of the following form. 

5 

For /i(S) = KeX^+y 7, where S is a dimensionless entropy, K is a constant, and h is a 
differentiable function of S 



Of — 1 ^ 

P 


—pu 


iT 

-pv p 


and the upper triangular part of the symmetric matrix Uw is 


(2.12b) 



apu apv 

apu 2 — p apuv 

apv 2 — p 


_i_£! 

7-1 p 


f p(u 2 +v 2 )-^p 
u[\p(u 2 + v 2 ) — bp] 
v[^p(u 2 + v 2 ) - bp] 

- bp{u 2 + v 2 ) + I p(u 2 + v 2 ) 2 


(2.12c) 


Here, p * and p are related through 


P 


x e = x(pp 7 ) ( “ +t) > 


(2.12d) 
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P = X *[(p*)W] lT H (2.12e) 

with x ~ The variable p* and /3 are given by 

= ( 2 . 121 ) 

P = (2-126) 

where a is a constant. The constants a, b and c are, a = (1 - a - 7)/a, b ~ 7/(7 - 1 ) and 

c = (1 - 2 7 )/(t - 1). 

The flux vectors, expressed in the W variables are given by 

F (V(W)) = ^+p* =ja -^( Wl + l E 2p.)] r , (2.12h) 

COW) = £[-»» £+P* -^(«’i + ^tp‘)] T . (2.12i) 

The upper triangular part of the symmetric matrices F(U(W)) w and G(U(W)) w , expressed in 
the 17 variables are given by 


*W= - 

P 


apu apu — p 
u(apu 2 — 3p) 


apuv 

v(apu 2 — p) 
u(apv 2 — p) 


w[fp(u 2 + v 2 )-6p] 

4- cpu 2 + f p(u 2 + v 2 )ii 2 - |p(u 2 4- r 2 ) 
uv[cp4- §p(u 2 + v 2 )] 
u[6c^ 4- cp(u 2 4- v 2 ) + fp(u 2 4- v 2 ) 2 ] 

(2.12j) 


apv 



apuv 

v(apu 2 — p) 


apv 2 — p 
u(apv 2 - p) 
v(apu 2 — 3p) 


r[fp(u 2 -f v 2 )-Jp] 

txv[cp4- §p(« 2 +v 2 )] 

— b 4- cpv 2 — |p(u 2 4- v 2 ) 4- § p(« 2 4- v 2 )t> 2 

v \bc^ 4- cp(u 2 + v 2 ) + ~p(u 2 + v 2 ) 2 ] 

(2.12k) 


In all of the numerical examples presented in Gerritsen (1996), a = 1 - 27. 


The high order base scheme using the SHOEC splitting applies the fourth and sixth-order 
central differencings to F m (and G y ) and W m . Note that this splitting of the flux consists of 
a conservative and a non-conservative part. The non-conservative part appears not to produce 
wrong shock locations traditionally associated with the use of non-conservative formulations of 
the Euler equation for computations. This splitting seems to require less numerical dissipation for 
the Euler computations over the non-split form. See Gerritsen (1996) for illustration. Recently, 
Vinokur (1998, private communication) extended the SHOEC idea to a thermally perfect gas and 
nonequilibrium flows (consisting of mixtures of thermally perfect gases). 
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Stable Boundary Schemes for High Order Base Schemes : Consistent and stable one-sided 
boundary stencils for numerical boundary treatments using fourth and sixth-order central spatial 
differencing have not been available until recently. These stable boundary schemes employ 
the summation-by-parts idea to derive an energy estimate for the high-order central spatial 
differencings as interior schemes. See Kreiss & Scherer (1977), Strand (1994), Gustafsson & 
Olsson (1995), Olsson (1995) and Gerritsen (1996) for the boundary scheme formulae. 

2.1.3. Spatial Discretizations for Viscous Terms 

For simplicity let V mm be a viscous term in one dimension. The possible high order base schemes 
for V zx can be 

Central Differencings : (fourth and sixth-order) 

V„. * (v i+2 - 16V j+1 + 30 V, - 16V}-! + j , (2.13) 

V” ~ 180 ^ 2 ( 2V >+* ~ 27V i+ 2 + 270V,+i - 490V, + 270V,_, - 27V,_ 2 + 2V,_,j . (2.14) 
Compact Central Differencingt: (fourth and sixth-order, Hirsh (1975), Ciment & Leventhal 


(1975) and Lele(1992)) 

V,.« -^(c.-'il.v) , (2.15a) 

where for a fourth-order approximation 

(C.V)i = i (v i+ , + lOVj + (2.15b) 

(D.V) } = V j+l - 2Vj + Vj-i , (2.15c) 

and for a sixth-order approximation 

{C x V)j = V J+1 + «o Vj -j- Vj~ i, (2.15d) 

(D.V), = bo (Vj +t - 2Vj + ^ (v j+2 - 2V, + Vj_ 2 ) , (2.15e) 

a 0 = 5.5, (2.15f) 

b 0 = 4(a 0 - l)/3, (2.15g) 

c 0 = (10 - a 0 )/3. (2.15h) 
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2.2. Characteristic Filters 

In this section we first discuss the procedure for applying the characteristic filter for multistage 
and LMM types of time discretizations. We then discuss forms of the characteristic filter. 


2.2.1. Procedure to Apply the Filter Step 


If a multistage time discretization like the Runge-Kutta method is desired, the spatial base 
scheme discussed in the previous section is applied at every stage of the Runge-Kutta method. If 
viscous terms are present, it is more consistent to use the same order and type of viscous base 
scheme as for the convection terms. After completion of the full Runge-Kutta step, a characteristic 
based filter step in the form of nonlinear dissipative numerical flux terms is applied at the end. 
Let U n+1 be the solution after the completion of one full Runge-Kutta step of the base scheme . 
advancement (i.e., replace U n+1 with (7 n+1 in (2.5) and (2.6)). Let Lf be the filter operator with 


L f {F\G')^ = 


Ax 


r»* p* 


1 

+ — 


g* . -G* 


Ajr p'* + » ~ *»-* 


(2.16) 


where F ?. , . and G* . , are the dissipative numerical fluxes for the filter operator. Hereafter, 
we refer to ‘ ‘F* x . ” and “GT . , , ” as the ‘‘filter numerical fluxes”. Possible formulae for the 

J i j f" Ji* ‘ j 

filter numerical flux will be discussed shortly. Then, we define the new time level as 


U n+ 1 = U n+1 + A tL f (F\ G*) jt k- (2.17) 

The filter numerical flux F* x . and G* , . t are evaluated at U n+1 . 

J+ ^ j 

If one desires a time discretization that belongs to the class of LMMs, then the filter operator 
L f can be applied as a dissipative numerical flux in conjunction with the base scheme. The filter 
numerical flux F* x , and G* , , in this case are evaluated at U n for explicit LMMs. For implicit 
LMMs additional similar filter numerical fluxes evaluated at the n + 1 time level are involved. 
Alternatively, procedure (2.17) can be applied to LMMs as well, where ?7 n+1 is the solution after 
the completion of one LMMs step of the base scheme. 

For time marching to steady states using implicit LMMs, certain flow physics only requires an 
explicit dissipation term. Also, the implicit operator can be different from the explicit operator. 
See Yee (1985), Yee et al. (1990) for some efficient conservative linearized implicit forms. 


2.2.2. The Filter Numerical Fluxes 

There are many possible candidates for the filter operator in conjunction with high order base 
schemes. Here, we propose to use filter operators that have similar width of grid stencils as the base 
scheme for efficiency and ease of numerical boundary treatment. Higher than third-order filter 
operators are of course applicable, but they are more CPU intensive and require special treatment 
near boundary points for stability and accuracy. On one hand, it would defeat the purpose of 
achieving efficiency. On the other hand, near shocks and shears, the resolution of higher than 



11 


third-order TVD or ENO schemes is comparable to their lower-order cousin except with a slight 
gain in resolution in regions near steep gradients and smooth flows. (Engquist & Sjogreen, 1995; 
Donat, 1994; Carpenter & Casper, 1997). If the lower-order filters are able to minimize but at 
the same time provide the proper amount of numerical dissipation away from shocks and shears 
to stabilize the non-dissipative nature of the high order base scheme, we would achieve similar 
resolution with improved efficiency. This is the philosophy used to choose the filter numerical 
fluxes. 

The simplest form is a scalar filter proposed by Gustafsson and Olsson (1995). It is the 
same form used by Jameson et al. (1981) to supply a second-order dissipation to a low order 
(second-order) central differencing for shock-capturing purposes. They used a switch similar to 
that of Harten (1978). The Harten switch was designed for self-adjusting hybrid schemes between 
Harten’s first-order ACM scheme and second or higher-order schemes. Instead of switching from 
a higher-order scheme to a first-order scheme for shock and shear capturing, we generalized 
Harten’s idea of achieving, in a loose sense, a low dissipative high order shock-capturing scheme 
by nearly maintaining the accuracy of the high order nondissipative property using characteristics 
based filters. The reason for the characteristic base filters is that scalar filters do not take into 
account the different wave characteristics of the Euler equations. For complex shock waves, shear 
and turbulence interactions, one has better control of the amount of dissipation associated with 
each wave. 


Filter Numerical Fluxes and Nonlinear Dissipation of Shock- Capturing Schemes : 

We start with any second or third-order TVD or ENO scheme that can be recast as the sum of 
central differencing and nonlinear dissipation terms. Recall that the Harten, Harten and Yee 
and Yee’s Symmetric TVD schemes (Harten, 1983, 1984; Harten & Yee, 1985; Yee, 1985) are 
already cast in this form. For example, let L tv d be a TVD (or ENO) spatial operator with 


1 

j _ ui . 

1 


Ltvd(F) G)j t h = 

Fj+ 

_|_ __ 

Ay 

G j,‘+§ ~ G i,*-i 


( 2 . 18 ) 


Take for example, the F flux. We cast the numerical flux F i+b h into the following form 

1 


F. 


,t = 2 




Fj+i,h + Fj,h + Rj+$ 

^i+f 


( 2 . 19 ) 


Here, | [Fj + i^ + Fj ,*] is the central differencing portion of the numerical flux F i+\+' and the 
last term R j+ i$ j+ i (with the suppression of the k index), is the nonlinear dissipation. For 
characteristic based methods, the quantity R j+ i is the right eigenvector matrix of ~ using, for 
example, the Roe’s approximate average state. Note that the eigenvector Rj + 1 should not be 

confused with the R in (2.2). Similarly, we cast the G jfh+ ± in the same manner. 

We use these nonlinear dissipation terms in Conjunction with Harten’s switch applied to each 
characteristic wave as the filter numerical fluxes. In essence, the nonlinear dissipation terms act 
as second or third-order ACM-like operators instead of Harten’s first-order ACM. The switch is 
used to signal the amount of nonlinear dissipation to be supplied to the high order nondissipative 
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scheme, one wave at a time. Thus, the current approach is also different in spirit from using ACM 
to sharpen the contact discontinuities in the original Harten second-order TVD scheme (Harten, 
1 983). Now we discuss the filter numerical flux in more detail. It is of the form 



2 R j+$*j+\' 


( 2 . 20 ) 


This filter numerical flux “F? +1 fc ” should not be confused with the standard numerical flux 

(( 2j ~ 2 * 

) ,h in (2.18). j is the modified form of the nonlinear dissipation portion of 

. The elements of denoted by 1 * are 

= "%+l+i+V ( 2 - 21 ) 

* n (2-21) are the elements of 4 in (2.19). The function 1 is the key mechanism for 
achieving high accuracy of the fine scale flow structure as well as shock waves in a stable manner. 
In other words, the elements of 1 are the same as the nonlinear dissipation term of the TVD or 

ENO scheme (2.19) with the exception of premultiplying by k6\ j. The parameter k is problem 
dependent. For the numerical examples in Section III, different examples require a different value 
of k. The range of k for these problems is 0.03 < k < 2. The function 0 l . + 1 is the Harten switch. 
For a general 2m -f- 1 points base scheme, Harten recommended 


3 = 






For all of the numerical examples, we use p — 1 and 


0 l . 


*+i 


max 



(2.22) 

u 

P 

(2.23) 

,1 1 


3+0 

!. 

(2.24) 


The i are elements of (tfj+i — Uj). The shock-turbulence interaction problems appear 

to favor this form of 0 l . . , . 

j+f 


Formulae for <j>\+ i are well known and can be found in the literature. For illustration purposes, 
we show a form of the * function in which all of the examples shown in Section III are used 
for the computations. We choose the Harten and Yee upwind TVD form where 



= ^(«y+ })(«}+, + 9j) - VK«' + 1 




(9j+i 

0 


«'+! # 0 


(2.25a) 

(2.25b) 
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Other shock-capturing schemes with structure similar to (2.19) such as the symmetric TVD 
schemes (Yee, 1985, 1987) and Roe second-order upwind scheme (Roe, 1983) are also applicable. 
The a* + 1 , l — 1, 2, 3, 4, are the characteristic speeds of evaluated at some symmetric average 
of Uj^h and Uj+ 1 ,*. The Roe’s average is an example of this (Roe, 1981). The function i/> is an 
entropy correction to |a*. x |. One possible form is (Harten & Hyman 1983) 

J-r 3 




+ Sl]/2Si 




(2.25c) 


For problems containing only unsteady shocks, is usually set to zero. Note that entropy-violating 
phenomena occur only for steady or nearly steady shocks. For steady-state problems containing 
strong shock waves, a proper control of the size of 6i is very important, especially for hypersonic 
blunt-body flows. See Yee et al. (1990) for a discussion. 

Examples of the commonly used ‘limiter’ function can be expressed as 


g) = minmod(a‘._ i,a'. + i), (2.25d) 

s‘i = / (“>+§ +a j-l)> (2.25e) 

9 l j = |o‘-i [(«' + ! )’ +^] +«' + i [(«'■_.) +^]}/[( a j+i) +( a ‘-i) + 2 *> 

9j = minmod(2a'._j,2a'. + ,, i(a' + | + <»'_.)), 
g\ - S - max 0,min(2|a^ + i |,5 • a^._i),min(|a^. + i | ? 25 • a^._i ) 


(-2.25f) 

(2.25s) 


; S = sgn(a'. , i). 


(2.25h) 


Here 6 2 is a small dimensionless parameter to prevent division by zero and sgn(a^. + 1 ) = 
sign(a I . i). In practical calculations 10" 7 < S 2 < 10 -5 is a commonly used range. For 
a 1 , t + — 0, g\ is set to zero in (2.25e). The minmod function of a list of arguments is 

equal to the smallest number in absolute value if the list of arguments is of the same sign, or is 
equal to zero if any arguments are of opposite sign. Later development in limiters have flooded the 
literature and has created much debate. Most of the improvements have been problem dependent. 
See Donat (1994), Engquist & Sjogreen (1995) and Jin & Lin (1996) on the error propagation 
for nonlinear approximations to hyperbolic equations containing discontinuities in derivatives or 
discontinuous solutions. 


2.3. Computer Implementation 

To avoid additional logical statements in the actual coding and to promote parallelization, 
several of the forms with the potential of dividing by zero are modified. They are: 


tj>(z) = + z 2 )’ 


(2.26a) 


14 


T'i+J 


! #»'•+. )(gj+i -gjK+i 

5 (“W 2 + e 


(2.26b) 


We use the switch 


,i = ll a i-M/2l l a i-i/2l 

l a jf+l /2 I + l a jf-l/2l + € 


(2.26c) 


In all of the computations, we take e = 10“ 7 . The value of S was taken to be 1/16 (unless 
indicated) to satisfy an entropy condition. However, the fine scale flow structure showed minor 
sensitivity to the value of this constant. 


2.4. Other Applicable Base Schemes 


There are other possible high order base schemes that one can use. For example, the fifth or 
seventh-order upwind schemes and the 2-4 or 2-6 MacCormack scheme. If the fifth or seventh- 
order upwind schemes are used as base schemes, one needs to subtract the dissipation portion 
of these upwind schemes from the filter step. The dissipation portion of these upwind schemes 
can be obtained by rewriting the scheme into two parts, a central part and the “rest”. The “rest” 
is the dissipation portion. For the 2-4 or 2-6 MacCormack scheme, the time discretization is an 
integral part and one has to use the complete scheme as the base scheme. For the filter step, one 
adds the filter numerical fluxes as an added corrector step as described in Yee (1989) or Yee & 
Shinn (1989). The filter numerical flux is the same as (2.20) but the <t> l j + i has a slightly different 

form to take into account the Lax-Wendroff type of At 2 term. See Yee (1989) or Young & Yee 
(1987) for the formula. 


2.5. Other Applicable Characteristic Filters 


MUSCL Approach Using an Approximate Riemann Solver: The filter numerical flux function 
for an upwind MUSCL-type scheme as described in Yee (1985a, 1989) using an approximate 
Riemann solver can be expressed as 



1 

2 


*»**i + r 


The elements of , and the vector (a‘) i+ i are given by 


(2.27a) 


(*°)' + . =-*(« 0 )' + ^(K)' +l )(“ 0 )y +§ . 



(2.27b) 

(2.27c) 


where ^>((a° )*. + 1 ) can be |(a° )*. + i 1 or the same form as (2.25c). Here (o° )*. + 1 are the eigenvalues 
and R°j + i is the eigenvector matrix of | ^ evaluated using a symmetric average between Uf+i 
and Uf, , ; i.e., 

3~y 2 
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(a°) l . + i =a\V* +v U]‘ +i ), (2.27d) 

^ + . =R(U« + ,,Uf + >). (2.27e) 

The switch (0°)* + i is the same as (2.22) and (2.23) except it is evaluated using a symmetric 

average between XJf i and Uf , . However, there are options in applying the limiters for system 

- ' — j-f- j 

cases. Namely, one can impose limiters on the conservative, primitive, or characteristic variables. 
The U R and U L are the upwind biased interpolation of the neighboring U^j, values with slope 
limiters imposed. 


MUSCL Approach Using the Lax- Friedrichs Numerical Flux : The filter numerical flux func- . 
tion for a MUSCL-type approach using the higher-order Lax-Friedrichs numerical flux (Yee, 1989, 
1997) can be expressed as 




(2.28a) 

where j is 





(2.28b) 

and (o° )™f can be 

(«TA' = *(K + )l+ c m). 

(2.28c) 


where | <X<1- The overbar for the quantity (U R + 1 — U r j* + i ) means the Harten switch together 
with k is applied to each element of the vector. However, the a *. + 1 in 8* + j are replaced with a 
jump in the the conservative variables U R t and Iff, »• If primitive or other variables are used 

3 3 < 3 

for the right and left states, the switch together with n should be applied to the corresponding 
variables. 


2.6. Filter Numerical Flux for Time-Marching to Steady States 

For time-marching to the steady states, one usually needs to add a fourth-order dissipation to 
a second-order spatial differencing scheme (Beam & Warming, 1978). For the present schemes 
using characteristic filters, in addition to the filter operator Lf, one might need to add a sixth-order 
dissipation to a fourth-order spatial base scheme and an eighth-order to a sixth-order spatial base 
scheme in regions away from shocks for stability and convergence. Let L * be such an operator. 
Take the case of a Runge-Kutta time discretization as in 2.2.1. There are two ways of incorporating 
the Ld operator. One way is to incorporate the L& operator at every stage of the Runge-Kutta 
method. Another way is to include the L operator as part of the filter step (2.17); i.e., 


u n+1 = U n+1 + A tLf(F',G')j, h + A tL d (V n+x )j <k , 


(2.29) 
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where Lf(F m )j t h is the same as in 2.17. These two ways of applying the Ld operator are most 
likely problem dependent and time integrator dependent. Extensive numerical experimentation is 
needed. For LMM type of time integrators, Ld is used in conjunction with Lf as an additional 
dissipation as discussed in Section 2.2. 1 . 

To minimize the amount of Ld in the vicinity of shock waves, there should be a switching 
mechanism (different from « in (2.21)) to turn off the Ld operator in the vicinity of shock 
waves. The Ld operator can be applied to the conservative, primitive or the characteristic variables. 
The simplest form is to apply Ld to the conservative variables. Alternatively, since all of the work 
in computing the average states and the characteristic variable are done for the Lf operator, one 
can apply the Ld operator on the characteristic variables instead of the conservative variable. In 
this case, the switching mechanism kd can be a vector so that it is more in tune with the shock 
detector of the approximate Riemann solver. 

III. Numerical Examples 

In all of the computations the classical fourth-order Runge-Kutta time discretization is employed. 
The detailed programming allows the Euler and viscous terms to be computed using separate 
methods. The basic spatial schemes are (i) non-compact central, (ii) compact central and (iii) 
predictor-corrector upwind or upwind biased. Non-compact schemes are the standard second, 
fourth and sixth-order methods. Compact schemes are either the standard symmetric fourth-order 
or the sixth-order Pade schemes. For the purposes of this paper we concentrate on the central 
schemes with the same order of accuracy and type of base scheme for the convection and viscous 
terms. Comparable accuracy was obtained with the upwind or upwind biased schemes proposed by 
Hixon and Turkel (1998). The filter operator (2.16) in conjunction with (2.20) - (2.26) are used as 
a filter step at the end of the full Runge-Kutta time step. Hereafter, we refer to this approach as the 
ACM/TVD (or simply ACM) method, indicating the fact that only one type of TVD dissipation 
is used for the numerical study. The various combinations of schemes considered for numerical 
experiments are shown on Table 3.1. The notation shown in Table 3.1 will be used for discussing 
the results for different numerical schemes. Here, the notation “TVD” with the various orders 
attached at the end means the second-order TVD dissipation (without the ACM switch) is used as 
the filter with the indicated order of the base scheme for the convection and viscous terms. For 
simplicity of discussion, unless otherwise indicated, the term TVD or ACM scheme means the 
selected base schemes indicated in Table 3.1 using the TVD or ACM/TVD filter. Studies using 
ENO dissipation as filters are planned. Computations using the SHOEC splitting in conjunction 
with high order central differencings as base schemes for a variety of perfect gas and nonequilbrium 
flow applications are also planned. It appears that the SHOEC splitting is more stable and requires 
less numerical dissipation. 

Without introducing additional notation, for inviscid flow simulations the same notation is used 
except the viscous terms are not activated. 

Five test cases are considered. The first two are inviscid and the last three are compressible 
DNS computations. These test cases are chosen to examine the versatility and accuracy of the 
proposed schemes for a variety of flows where most conventional methods exhibit difficulty in 
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obtaining low dissipative solutions in an efficient and stable manner. All the test cases use either a 
uniform or mildly stretched Cartesian grid in the ^-direction. No attempt is made to enhance the 
resolution using appropriate adaptive grid procedures. At present, the code used for the test cases 
reduces to lower order central base schemes near the boundary points. Stable boundary treatments 
suggested by Gustafsson and Olsson (1995) should be used and are not yet implemented for the 
present study. Nonreflecting boundary conditions or characteristic inflow and outflow boundary 
treatment are also not implemented. 

The five test cases are: ( 1 ) a stationary vortex evolution, (2) a horizontally convecting vortex, (3) 
a vortex pairing in a time-developing mixing layer with shock waves formed around the vortices, 
(4) a shock wave impinging on a spatially evolving mixing layer where the evolving vortices 
must pass through a shock wave, which in turn is deformed by the vortex passage, and (5) a 3-D 
compressible turbulent channel flow to validate that the proposed schemes are in fact capable of 
sustaining turbulence. To examine the resolution of the proposed schemes where shock waves are 
absent, the computation was compared with the CEN44 (the classical spatially fourth-order central 
differencing for the convection and diffusion terms) before shock waves were developed for the 
vortex pairing case. Good agreement was obtained. 

Aside from evaluating the vortex preservation property, the performance of these schemes with 
the presence of shock waves and turbulence are evaluated based on the following factors: 

(a) Effect of the ACM term 

(b) Effect of the order of the base scheme 

(c) Effect of the grid size (grid refinement study) 

(d) Effect of employing a compact or non-compact base scheme 

(e) Effect of the adjustable parameter k for the particular physics 

(f) Effect of the flux limiters 

(g) Shear and fine scale flow structure capturing capability 

3.1. Isentropic Vortex Evolution 

The first two test cases are chosen to assess the performance of the proposed schemes for 
evolution of a 2-D inviscid isentropic vortex in a free stream. Similar test cases have also been 
used by several authors for testing other schemes with respect to vortex preservation (Gerritsen, 
1996; Davoudzadeh et al., 1995; Shu, 1996). The mean flow velocity, «oo and t>oo, pressure, poo, 
and density, p are considered to be free stream. Test case 1 is a stationary (steady) vortex with 
(woo,Uoo) = (0,0), and case 2 is a horizontally convecting vortex with (uoo?Uoo) = (1,0). In both 
cases poo = poo = 1. 

As an initial condition, an isentropic vortex with no perturbation in entropy (SS = 0) is added 
to the mean flow field. The perturbation values are given by 

(Su,Sv)= (~y,x), (3.1) 


(7-I)/? 2 . r . 
SjTT 2 ’ 


6T = 


(3.2) 
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where ft is the vortex strength and 7 = 1.4. Note that the vortex strength ( 3 should not be confused 
with the ft in section 2.1.2. Here T = Too = 1.0, (*,y) — (x — *„ 0 ,y — y Vo ), where a? Vo and 
y„ 0 are coordinates of the center of vortex initially, and r 2 = x 2 -f y 2 . The entire flow field is 
required to be isentropic so, for a perfect gas, — 1. 

From the relations, p = /><» + 8p, u = -f 8u, v = v^, + 8 v, T — Too + £T, and the above 
isentropic relation, the resulting state for conservative variables is given by 


(3.3) 

(3.4) 

(3.5) 

(3.6) 

(3.7) 

These two vortex problems provide a good test bed for evaluating the schemes performance 
with the absence of shock waves and turbulence. The exact solution with given initial states is just 
a passive convection of the vortex with the mean velocity and thus provides a good measure of the 
accuracy of the schemes for smooth solutions of the Euler equations. Figure 3.1 shows the initial 
vortex covering a domain of 0 < x < 10 and — 5 < y < 5. 

Periodic boundary conditions (BCs) in both directions are traditionally used for these test 
cases. Since the code reduces to lower-order central base schemes near the boundary points, and 
nonreflecting BCs are not used, non-periodic BCs simulations would provide the opportunity to 
examine the effect of the sizes of computational boxes on preserving the vortex core during time 
evolution. Nonreflecting BCs and stable boundary treatments suggested by Gustafsson and Olsson 
(1995) will be implemented for a future study. 

Both test cases employ a uniform Cartesian grid. Density profiles at the centerline, y — 0, 
cutting through the center of the vortex of the various schemes are used for comparison. Data on 
the centerline was extracted up to 5 unit lengths away to the left and the right, from the location 
of the center of the vortex. In all of the computations for the vortex evolution, unless otherwise 
indicated, 8 = 0.01 (2.264), and limiter (2.25f) are used. 

3.1.1. Stationary Vortex 

For the stationary vortex test case, a uniform grid spacing of Ax = Ay = 0.125, covering the 
domain of 0 < x < 50 and -5 < y < 5 are used. The grid is 401 x 81. The vortex is placed at 
the center of the rectangle, (x - x V(n y - y Vo ) = (25,0). For reasons of economy, only the left 
and right boundaries in the aj-direction are kept fairly distant from the center of the vortex core at 
25 unit lengths. Only 5 unit lengths are used in the y-direction. Since there are no shock waves or 


p = (Too + 8T)t-' = 
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steep gradient regions for this flow, the filter is used to stabilize the nonlinear governing equations. 
For this reason, the filter coefficient k (2.21) should be kept very small. We use k = 0.05 for all 
of the computations. The time step is also fixed with At — 0.04. If we use the SHOEC method as 
the base scheme suggested by Gerritsen (1996), the filter step is not necessary (see Section 2.1.2 
for a discussion). Studies using the SHOEC method (Gerritsen, 1996) as the base scheme and the 
accompanying stable boundary treatment of Gustafsson & Olsson (1995) are planned. 

The density profiles across the vortex at the centerline, y = 0, for the various schemes are 
illustrated in Fig. 3.2. Figure 3.2 which shows the effect of increasing the accuracy from second 
to fourth and sixth-order using the TVD filter (TVD22, TVD44, and TVD66) compared with the 
ACM/TVD filter (ACM22, ACM44 and ACM66). Although the order of the viscous terms are 
indicated on the method, the viscous terms are not activated. Figure 3.2a compares the exact 
solution with the solutions obtained by the TVD22, TVD44 and TVD66 methods at t — 50 (after 
1250 times steps). All of the TVD methods, regardless of the order of the base schemes are very 
diffusive, especially around the vortex core. Higher-order base schemes exhibit slightly better 
resolution than the second-order base method. Figure 3.2b displays the same computation at a later 
time t — 100 (after 2500 time-steps). The computed vortex core is even more diffused compared 
with the exact solution. 

Figure 3.2c,d shows the same computation at t = 50 and 100 using the various orders of the 
base scheme with the ACM/TVD filter. These figures display the effect of the ACM/TVD filter 
on the vortex core. The ACM methods, regardless of the order, have not diffused the vortex core 
after t — 50. All numerical solutions fall almost on top of the exact solution, except for very small 
differences for the ACM66 method. At t ~ 100, the ACM66 resolution has been slightly displaced 
due to the boundary effects. However, the ACM22 and ACM44 remain quite accurate. 

3.1.2. Horizontally Convecting Vortex 

For the horizontally convecting vortex, again, a uniform grid spacing of A* = Ay = 0.125, 
covering the domain of 0 < * < 110 and -5 < y < 5 is used. The grid is 481 x 81 for t = 50 
and 881 x 81 for t = 100. The vortex is initially placed at (* - x Vo1 y ~ y Vo ) = (5,0). The time 
step, At = 0.04, is fixed for all runs, as is the vortex strength, 0 = 5. The adjustable parameter, 
k, is set equal to 0.05 as in case 1 . The vortex is convected to the right by the mean flow velocity. 

The physics of the present vortex evolution is similar to the stationary case, except the vortex 
is convecting. Since the ACM44 and ACM66 are less diffusive than the ACM22 for case 1 , only 
the ACM44 and ACM66 are used for the present computation. Figure 3.3 displays density profiles 
across the centerline at y = 0, comparing the exact solution with the ACM44 and ACM66 methods 
at t = 50 and t — 100, and after 1250 and 2500 time steps respectively. All numerical solutions 
are very accurate and fall almost on top of the exact solution. In these computations no visible 
boundary effects are seen because the right boundary of the domain in the direction of the vortex 
convection is initially kept relatively far away. Although not visible from the density profiles, 
ACM44 exhibits small oscillatory solutions at t — 100. However, the ACM66 exhibits only small 
oscillations at the outer edge of the vortex. Figure 3.4 shows the density contours comparison of 
the exact solution with ACM44 and ACM66 at t - 100. We reran the same case using ACM44 and 
an increased k — 0.07. The small oscillation disappeared and the solution is as accurate as for the 
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ACM66 using k = 0.05. Both ACM44 and ACM66 exhibit good vortex preservation. Figure 3.4 
also shows the comparison of the two different values of k using ACM44. These results indicate 
the effect of k on the various orders of the ACM methods. For viscous flow, in the presence of 
shocks, shears and turbulence, the effect of k on the resolution of fine scale flow structure plays a 
different role than for the in viscid flows with smooth solution. When we ran the same case using 
CEN44, the solution blew up at t — 65. 

3.1.3. Boundary Effects on the Stationary Vortex 

The effect of the size of the computational boxes was studied for the stationary vortex evolution, 
case 1 . All the numerical experiments with the TVD and ACM methods discussed for case 1 were 
repeated on a smaller computational domain of 0 < * < 10 and — 5 < y < 5 for which the initial 
vortex touches the boundary of this smaller computational box. The grid spacings and time-step 
are the same as before, except the grid is now 81 x 81. Figure 3.5 shows the computations on 
the reduced domain. Comparison with Fig. 3.2 clearly demonstrates the effect of the physical 
boundary distances in the ^-direction. The discrepancies between corresponding results on the 
larger and smaller domains are more pronounced at t = 100 for the ACM methods. Figure 3.6 
compares the numerical solution of the ACM66 method on larger and smaller domains with the 
exact solution. 

Figure 3.7 displays the effect of the adjustable parameter k in controlling the boundary effects 
(aj-direction) for the ACM44 and ACM66 methods at t- 100. We reduced the k value from 0.05 
to 0.035 on the same smaller computational box. The profile for both methods are on top of each 
other. The deviation from the exact solution of the computed solution due to boundary effects is 
less pronounced than for k = 0.05. Figure 3.8 shows the effect of using different limiters (limiter 
(2.25f) vs. limiter (2.25h)) at* = 100 on the same smaller computational box. 


3.2. Vortex Pairing in a Time-Developing Mixing Layer 

This test case studied vortex growth and pairing in a temporal mixing layer at a convective 
Mach number equal to 0.8. At this Mach number there are shock waves (shocklets) that form 
around the vortices and the problem is to compute accurately the vortex evolution while avoiding 
oscillations around the shocks. Previous calculations of the problem can be found in Sandham 
and Reynolds (1989), Lumpp (1996a,b) and Fu and Ma (1997). Here we set up a base flow as in 
Sandham and Yee ( 1 989) 

u = 0.5tanh(2j/), (3*8) 

with velocities normalized by the velocity jump u\ — «2 across the shear layer and distances 
normalized by vorticity thickness. 


6 


w 


U l ~ U 2 
(du/dy) 

max 


(3.9) 


Subscripts 1 and 2 refer to the upper (y > 0) and lower (y < 0) streams of fluid respectively. The 
normalized temperature and hence local sound speed squared is determined from an assumption of 
constant stagnation enthalpy 
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C 2 = c\ + Iyi(«; - «*). (3.10) 

Equal pressure through the mixing layer is assumed. Therefore, for this configuration of u 2 — — u\ 
both fluid streams have the same density and temperature for y — >• ±oo . The Reynolds number 
defined by the velocity jump, vorticity thickness and kinematic viscosity at the free-stream 
temperature is set here to be 1000. The Prandtl number is set to 0.72, the ratio of specific heats 
is taken as 7 = 1.4 and Sutherland’s law with reference temperature Tr = 300JfiT is used for 
the viscosity variation with temperature. The reference sound speed squared, c] v is taken as the 
average of c 2 over the two free streams. 

Disturbances are added to the velocity components in the form of simple waves. For the normal 
component of velocity we have the perturbation 

2 

v = ah cos(2 irkz/L m + <£*)exp(-y 2 /&), (3.11) 

fc=i 

where L x = 30 is the box length in the ^-direction and b = 10 is the y-modulation. In our test 
case we simulate pairing in the center of the computational box, by choosing the initially most 
unstable wave Jfe — 2 to have amplitude a 2 = 0.05 and phase <j> 2 = —v/2, and the subharmonic 
wave k — 1 with ai = 0.01 and <j> 1 = — ir/2. The u — velocity perturbations are found by assuming 
that the total perturbation is divergence free. Even though these fluctuations correspond only 
approximately to eigenfunctions of the linear stability problem for a compressible mixing layer, 
they serve the purpose of initiating the instability of the mixing layer and have the advantage as a 
test case in that they can be easily coded. 

Numerically the grid is equally spaced and periodic in the *— direction and stretched in the 
y— direction, using the mapping 

L y sinhf&yT?) 

y= 

where we take the box size in the y-direction L y = 100, and the stretching factor b y 
mapped coordinate r] is equally spaced and runs from —1 to +1. The boundaries at 
taken to be slip walls. For example, at the lower boundary 


pi — P 2 i (3.13a) 

(/**)i = (/"Oil (3.13b) 

(/>t>)i = 0, (3.13c) 

(e). = [4(e) 2 - («),]/3, (3.13d) 


(3.12) 

= 3.4. The 
±L y / 2 are 


where subscripts here refer to the grid point and e is the total energy. 
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3.2.1. Computational Results 

We compute this test case on 401 x 401, 201 x 201, 101 x 101 and 41 x 41 grids. There 
is little to choose in the shock resolution properties with the variation in order of accuracy of 
the scheme since the proposed method will not improve the resolution of the shock waves over 
standard TVD or ENO schemes. We choose to compare temperature contours, which are most 
sensitive to oscillations (Lumpp, private communication), and accuracy of the fine scale structure. 
Figure 3.9 shows a snapshot of the temperature contours at t — 40, 80, 120 and 160 using ACM44, 
illustrating the roll-up of the primary vortices followed by vortex merging. Shock waves and 
shears form around the vortices with a peak Mach number ahead of the vortex of approximately 
1.55 at t = 120. The grid is 201 x 201. The majority of the comparisons, however, use a 101 x 101 
grid. In all of the computations for the vortex pairing case, unless indicated, limiter (2.25f), and 
6 = 1/16 (2.26a) are used. The majority of the computations used « = 0.7 (2.21) for the nonlinear 
fields for the ACM methods. 

It is noted that a similar vortex pairing was used by Shu et al. ( 1 99 1 ) to evaluate the performance 
of high order ENO schemes. The present results show superior performance over the result of Shu 
et al. 

Effect of the ACM Term and the Order of the Base Scheme : 

Figure 3.10 shows the effect of increasing the accuracy from second to fourth and sixth-order using 
the TVD filter (TVD22, TVD44 and TVD66). As can be seen there is almost no improvement 
as the order of accuracy is raised. Figure 3.1 1 shows the same plot using the ACM/TVD filter 
(ACM22, ACM44 and ACM66). Here there is an improvement, although the results even for 
the lowest order are quite good. All the ACM schemes capture the shock waves with minimal 
oscillations. Although not shown, the temperature contours for the TVD filter of various orders 
using a 101 x 101 grid are not nearly as accurate as the ACM44 using a 41 x 41 grid. See the 
last plot of Fig. 3.17 or Sandham & Yee (1998) for illustrations. It can be seen that there is a 
significant advantage in moving from second to fourth order, but there is a smaller gain in moving 
from fourth to sixth order using TVD or ACM/TVD as a filter. (This is contrary to the isentropic 
vortex convection where there are definite benefits of moving from fourth order to sixth order. The 
effect of order of accuracy are more pronounced for long time integration of pure convection.) The 
fine scale flow structures are nearly resolved using a 101 x 101 grid compared with the reference 
solution using ACM44 and a 201 x 201 grid. Results from the ACM method are far superior 
to those from the standard TVD formulation. Note that there is no improvement in the shock 
resolution among the various orders of the base schemes, since the ACM term limited the amount 
of dissipation away from shocks and steep gradient areas whereas the shock resolution is dictated 
by the flux limiter. 

Effect of the Grid Size (Grid Refinement Study): 

To investigate the effect of order of the accuracy in more detail we consider simulations on a 
very coarse grid of 41 x 41 points. Such a case corresponds in practice to simulation of scales 
of turbulence arising from shear layers only two or three computational cells across. Figure 3.12 
shows results for the ACM44 schemes. To ensure that the fine scale flow structure is fully resolved 
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by the 201 x 201 reference grid, the same simulation was done on a 401 x 401 grid (figures not 
shown). The resolution of the ACM methods on a 41 x 41 grid is comparable to TVD methods on 
a 101 x 101 grid. 

Effect of Compact or Non-Compact Base Scheme : 

For wave propagation and computational problems the performance of fourth and sixth-order 
compact schemes, although more CPU intensive, appears to be superior over their non-compact 
cousin. For problems with shock waves the benefit of compact over non-compact schemes is less 
known due to the filter step. For this purpose we conduct similar experiments using (2.9) as base 
schemes. Results for the fourth and sixth-order compact schemes are similar to results from the 
sixth-order non-compact scheme. Again there is little improvement compared with the fourth-order 
non-compact scheme. The compact schemes are almost double the CPU over their non-compact 
cousin for this 2-D compressible DNS computation using the fourth-order Runge-Kutta method. 
For this particular flow physics, a conclusion is that the use of the ACM in the filter step is essential 
to get the benefits of moving from second to fourth order, but even with the compact method as 
the base scheme, there is little benefit in moving to even higher-order base schemes. See Sandham 
& Yee (1998) for illustrations. This is in contrast to the isentropic vortex convection where there 
are benefits of moving from fourth order to sixth order for long time integration. 

Effect of the Adjustable parameter n and Shear and Fine Flow Structure Capturing Capability : 

The ACM switch has been demonstrated to give good shock resolution and to be essential if the 
benefits of higher-order discretization schemes are to be realized. There is, however, an adjustable 
parameter k in the formulation, and results are sensitive to the precise choice of its value. Figure 
3.13a and Fig. 3.13d illustrate the effect on the result using ACM44 for the pairing test case of 
reducing the parameter from 0.7 to 0.35. The vorticity and momentum thickness development is 
improved due to the reduction in numerical dissipation. From the temperature contours on Fig. 
3.13a,d it can be seen that this has been achieved at the cost of formation of small oscillations 
around the shock wave. For the present problem one would be ready to pay this price to get the 
more accurate vortex evolution. However, in general it is not known how such numerically-induced 
oscillations interact with small scales of turbulence. For the current method the correct procedure 
for a simulation of shock-turbulence interaction would be to find the smallest value of k to 
resolve the shock waves satisfactorily and then increase the grid resolution until the turbulence is 
adequately resolved. There are perhaps other formulations of the ACM switch parameter k based 
on the flow physics that can perform the adjustment to higher values automatically when stronger 
shock waves are present. This is a subject of future research. 

To balance the shear and shock-capturing, one alternative is to switch to a more compressive 
limiter (see Yee (1989)) for the linear characteristics fields. Another alternative is to reduce the 
value of k for the u,v linear fields. The comparison of using different values of n for the linear 
and nonlinear fields is also shown in Figs. 3.12b,c,d and Fig. 3.13b. 

Shock-capturing schemes are designed to accurately capture shock waves, but with a less 
accurate capturing capability for contact discontinuities. In fact, the mixing layer seen at a large 
scale is a contact discontinuity. If one uses enough grid points to resolve the region of high shear in 
conjunction with the physical viscosity, it might not need to be ‘captured’. Contact discontinuities 
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relate to the characteristic velocities u and v. As an experiment, these linear characteristics were 
computed with no numerical dissipation. Filters were only applied to the nonlinear characteristic 
fields u ± c and v ±c using the ACM44 method. Interestingly, the computation was no less stable 
than that of the full TVD or ACM schemes. Figure 3.13c shows the resulting temperature field. 
It can be seen that good results were obtained, although there is a trace of oscillation near the 
shock wave. This could be remedied by increasing « slightly. However, the flow features of the 
shear and fine flow structure are accurately captured with similar resolution as the 201 x 201 grid 
with ACM44 applied to all of the characteristic fields. Figure 3.14 shows the comparison of the 
different n values for the nonlinear fields with the « = 0.0 for the linear fields. 

Effect of the different flux limiters and the Adjustable Constant 8 (2.26a): 

Figure 3.15 shows the comparison of the five classical flux limiters (2.25d-2.25h) (see Fig. 
3.13a for lim3). Lim4 (2.25g) appears to perform the best with the A shock, shear layer and 
fine scale structure similar to the reference solution but using half of the grid size in each spatial 
direction. Figure 3.16 shows the effect of the various 8 ( (smu) 2 indicated on the plots) values 
on the fine scale flow structure capturing. It appears that 8 = 0 perform the best except in this 
case the TVD filter is technically entropy violating. For this computation « is set to 0.7 for all 
characteristic fields. 


3.2.2. First-Order Upwind Dissipation as Characteristic Filter Computations 


To examine the performance of using first-order dissipation as characteristic filters for high 
order base schemes, we considered an entropy modification of Roe’s first-order dissipation by 
redefining <f > 1 . , i in (2.21) to be 

J'r 3 





i+i 


(3.14) 


Figure 3.17 shows the computation using the modified Roe’s first-order dissipation as a filter with 
the various orders of central differencing base scheme denoted by ACM22/lst, ACM44/lst and 
ACM66/ 1 st. The 6 is set equal to 1/16 using (2.26a) for The result is comparable to 

TVD22, TVD44 and TVD66. Rerunning the same case with a smaller 8 results in more accuracy 
than using 8 — 1/16. 


3.2.3. Lax- Wendroff Type TVD Dissipation as Characteristic Filter Computations 


To examine the applicability of using the Lax-Wendroff type TVD dissipation filter, we 
recompute Fig. 3.1 1 using 
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In other words, using the Lax-Wendroff TVD dissipation as a filter involves an extra ^(a *. + 1 ) 2 

term in the t/> function compared with (2.25). This term is due to the U t t term in the Lax-Wendroff 
formulation. Note that the use of (3.15) is not consistent in two ways. First, the base scheme uses 
Runge-Kutta time integrator and the filter step uses the Lax-Wendroff type of mixed space and 
time formulation. Second, for 2-D, one has to use the Strang (1968) type of splitting to retain 
second-order time accuracy. 

Disregarding the inconsistency, Figure 3.18 shows a comparison of the various orders of base 
schemes using (3.15) together with the ACM switch, denoted by ACM22/i/ 2 , ACM44/t/ 2 , and 
ACM66A/ 2 . The solutions are more diffusive but slightly more stable than the corresponding cases 
without the v 2 term (Fig. 3.1 1). Although the use of (3.15) is not consistent, the results are far 
superior than the TVD22, TVD44 and TVD66 method (Fig. 3.10). 

3.3. Shock Wave Impingement on a Spatially-Evolving Mixing Layer 

The fourth test case has been developed to test the behavior of the schemes for shock waves 
interacting with shear layers where the vortices arising from shear layer instability are forced to 
pass through a shock wave. An oblique shock is made to impact on a spatially-developing mixing 
layer at an initial convective Mach number of 0.6. The shear layer vortices pass through the 
shock system and later through another shock, imposed by reflection from a (slip) wall at the 
lower boundary. The problem has been arranged with the Mach number at the outflow boundary 
everywhere supersonic so that no explicit outflow boundary conditions are required. This allows 
us to focus on properties of the numerical schemes rather than on the boundary treatment. 

Figure 3.19 illustrates the nature of the flow on a 641 x 161 grid illustrating the pressure, density 
and temperature fields using the ACM44 method and limiter (2.25f) with k = 0.35 for nonlinear 
characteristic fields and k = 0.175 for the linear characteristic fields. The parameter 6 is set to 
0.25. This computation is used as the reference solution. An oblique shock originates from the top 
left hand corner and this impacts on the shear layer at around x — 90. The shear layer is deflected 
by the interaction. Afterwards we have a shock wave below the shear layer and an expansion fan 
above it. The shock wave reflects from the lower solid wall and passes back through the shear 
layer. The lower wall uses a slip condition so no viscous boundary layer forms and we focus on the 
shock-wave interaction with the unstable shear layer. The full no-slip problem would, however, 
make a challenging test case for the future. 

The inflow is specified again with a hyperbolic tangent profile, this time as 

u = 2.5 + 0.5 tanh(2y), (3.16) 

giving a mixing layer with upper velocity = 3, lower velocity u 2 = 2, and hence a velocity 
ratio of 1.5. Equal pressures and stagnation enthalpies are assumed for the two streams, with 
convective Mach number, defined by 


M c = (3.17) 

Cl -I- c 2 

where ci and c 2 are the free stream sound speeds equal to 0.6. The reference density is taken as the 
average of the two free streams and a reference pressure as (pi -1- /> 2 )(^i — u 2 ) 2 / 2. This allows 
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one to compute the inflow parameters as given in the first two columns of Table 3 . 2 . Inflow sound 
speed squared is found from the relation for constant stagnation enthalpy ( 3 . 10 ). 

The upper boundary condition, given in column 3 of Table 3 . 2 , is taken from the flow properties 
behind an oblique shock with angle /3 = 12°. The table also gives the properties behind the 
expansion fan (column 4) and after the oblique shock on the lower stream of fluid (column 5) 
computed by standard gasdynamics methods with (3 = 38.118°. In practice, the conditions in 
regions 4 and 5 do not correspond exactly to the simulations due to the finite thickness of the 
shear layer. The Mach number of the lower stream after this shock is approximately M 5 = 1.6335 
and remains supersonic through all the successive shocks and expansion fans up to the outflow 
boundary. The resulting shock waves are not strong, but tests showed that they could not be 
computed without using shock-capturing techniques. The lower boundary was specified with the 
same slip condition used for the pairing case (Equation ( 3 . 13 )). 

The Prandtl number and ratio of specific heats were taken to be the same as for the vortex pairing 
test case. The Reynolds number was chosen to be 500. 

Fluctuations are added to the inflow as 

2 

v' = ^ ah cos(2nkt/T + fo)exp(-y 2 /6), (3.18) 

fc= l 

with period T = \/u c , wavelength A = 30, convective velocity u c = 2.68 (defined by 
u c = (ti 1 c 2 +« 2 Ci)/(ci +c 2 )) and b = 10. For k — 1 we take aj = 0.05 and<f> — 0, and for k = 2 
we take a 2 = 0.05 and ^ = ir/2. No perturbations are added to the u— component of velocity. 

The grid is taken to be uniform in * and stretched in y according to equation ( 3 . 12 ) with b y — 1. 
This stretching is much milder than for the pairing problem, as we have to resolve the shear layer 
even when it deflects away from y = 0. The box lengths were taken to be L 9 = 200 and L y — 40. 

The reference solution indicates that vortex cores are located by low pressure regions and the 
stagnation zones between vortices by high pressure regions. The shock waves are seen to be 
deformed by the passage of the vortices. Another interesting observation is the way the core of 
the vortex at x — 148 has been split into two by its passage through the reflected shock wave. 
In spite of the relatively high amplitude chosen for the subharmonic inflow perturbation there is 
no pairing of vortices within the computational box. We do, however, begin to see eddy shock 
waves around the vortices near the end of the computational box where the local convective Mach 
number has increased to around 0.66. The oscillations seen near the upper boundary for x > 120 
occur where the small Mach waves from the initial perturbations arrive at the upper boundary. The 
use of characteristic boundary conditions should remove this problem. Practically, the amplitude 
of oscillations is not sufficient to cause numerical instability or affect the remainder of the flow. 

The test case was also run on a grid of 321 x 81 with At = 0.12, 6 = 1/16, and k = 0.35 up to 
t = 120. Figure 3.20 and 3.21 show the effect of increasing the accuracy of the base scheme using 
the TVD and ACM/TVD filters. Again it can be seen that the ACM switch is essential for obtaining 
good vortex evolution (additionally better shock resolution is obtained). For a more quantitative 
comparison, see Figs. 18 and 19 of Sandham & Yee ( 1998 ). From Fig. 3.20 it is apparent that all 
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the standard TVD schemes, of whatever order, miss the correct vortex formation. From Figs. 3.20 
and 3.21 there is a visible benefit in moving from second to higher-order differencing, both in the 
amplitude of the fluctuations and in the correct convective velocity of the vortices. Figure 3.22 
shows the effect of the 6 values on the fine scale flow structure capturing. For this computation 
k = 0.35 is used for all characteristic fields. 

Disregarding the inconsistency. Fig. 3.23 shows a comparison of TVD44/i/ 2 with ACM44/i/ 2 on 
two grids. The solutions are again more diffusive but slightly more stable than the corresponding 
cases without the i/ 2 term. Although the use of (3.15) together with the ACM switch is not 
consistent, the results are far superior to TVD44 using (3.15) or (2.25) without the ACM switch. 

3.4. Compressible Turbulent Channel Flow 

As a final topic we consider the numerical simulation of three-dimensional turbulence. Given 
the additional numerical dissipation of shock-capturing schemes, even where shock waves aren’t 
present, it seems prudent to use this test case to validate that the new schemes proposed here are in 
fact capable of sustaining turbulence, before moving to complete simulations of shock-turbulence 
interactions. The case that we choose is similar to that of Coleman et al. (1995). In this problem we 
choose to normalize distances with half the separation between the walls, densities with the average 
density p, velocities with the friction velocity u T and temperature with the fixed wall temperature 
T w . We take the pure pressure-driven flow, where forcing terms corresponding to dp/dx — — 1 
are added to the a-momentum equation. For a computational domain we use a box size with a 
streamwise length of 3 (channel half widths), and spanwise length of 1.5. While it is not sufficient 
for the two-point correlation to fall to zero, these box lengths are adequate for common turbulence 
statistics, as illustrated on Table 3.3 which compares the standard incompressible results of Kim et 
al. (1987) with incompressible data computed by the first author with a fully spectral code on the 
small grid. It is likely that the same size domain will be a suitable test case for compressible flow 
provided the Mach number is low. At higher Mach numbers it is known (Coleman et al., 1995) 
that correlation lengths increase. The advantage of the small domain calculation as a test case for 
numerical methods is the small memory requirement making it feasible to run the calculations on 
workstations with limited memory. 

For the compressible calculation we take the Mach number based on the friction velocity and 
sound speed corresponding to wall temperature to be 0.05, giving a centerline Mach number of 
approximately 1 .1 . A full grid-refinement study has not yet been attempted for this configuration. 
Results have been prepared for a grid of 32 x 81 x 32 which is the same as that which has proved 
sufficient for good statistics up to triple moment budgets for the fully spectral code (Sandham & 
Howard, 1997). The ACM44 method was applied with « = 0.35 as in previous sections. The 
turbulence was indeed able to sustain itself, allowing the accumulation of statistical quantities for 
the flow. Figure 3.24a shows the mean velocity profile compared with the standard law of the wall 
(with Karman constant 0.41 and additive constant 5.5). The shift of the curve upwards relative 
to the incompressible result is consistent with (Coleman et al., 1995). Figure 3.24b shows the 
Favre-averaged stress -pu"v" and the total stress found by adding the contribution pdu/dy. The 
total stress must equal the non-dimensional y coordinate for the statistically converged flow. Favre 
averages are defined using mass weighting as u = pu/p and fluctuations given by u" = u - u. 
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3.5. Computational Costs 

For the DNS computations using the fourth-order Runge-Kutta method, the non-compact base 
schemes with the ACM/TVD filter are only around 25% more expensive than the same base 
schemes without ACM/TVD filter. This has been achieved by only requiring one application of 
the ACM/TVD filter per full time step for the convection terms. For LMM time discretizations, the 
non-compact base scheme with the ACM/TVD filter is only 10% more expensive than standard 
second-order TVD schemes. 

For the Cray C90 it was found that the compact schemes were significantly more expensive, but 
this result was distorted for the present code by incomplete vectorization. An extra cost of around 
2/3 is expected from considerations of operations count. 

Conclusions 

A generalization of the work of Gustafsson and Olsson and the ACM switch of Harten to a 
class of versatile low dissipative high order shock-capturing methods using characteristic based 
filters is proposed. The design principle of these schemes consists of two steps. The first step is 
the high order spatial and temporal base scheme. A variety of standard high order non-dissipative 
or low dissipative base schemes fits in the present framework. The second step is the appropriate 
filter for stability and shocks, shear-layers and fine scale flow structure capturing. Many of the 
TVD and ENO dissipations, after a minor modification, are suitable candidates as filters. The final 
grid stencil of these schemes is five if second-order TVD schemes are used as filters and seven if 
second-order ENO schemes are used as filters for a fourth-order base scheme. Numerical boundary 
condition treatment is simple and can be the same as for the existing base and filter schemes. 

The reason for proposing filter operators that have similar width of grid stencils as the base 
scheme is for efficiency and ease of numerical boundary treatment. Higher than third-order filter 
operators are applicable, but they are more CPU intensive and require special treatment near 
boundary points for stability and accuracy. It is well-known that near shocks and shears, the 
resolution of higher than third-order TVD or ENO schemes is comparable to their lower-order 
cousin except with a slight gain in resolution near steep gradients and smooth flows. The five 
test cases show that an increased in resolution with improved efficiency can be accomplished if 
we limit the proper amount of numerical dissipation away from shocks and shears to stabilize the 
non-dissipative nature of the high order base scheme. 

The approach is prompted partly by a need for an efficient method which is capable of 
highly resolved DNS of compressible turbulence in the presence of shock waves for a variety 
of flow speeds and partly by the need for the preservation of vortex convection and fine scale 
flow structure capturing. The five test cases illustrate the versatility of the proposed schemes in 
accurately capturing a variety of flow features, where most conventional methods exhibit difficulty 
in obtaining low dissipative solutions in an efficient and stable manner. Higher accuracy is 
achieved with fewer grid points when compared with standard TVD or ENO schemes. In all of 
the test cases, if adaptive grid refinements were used, additional gain in efficiency and accuracy 
should be realized. Application of these schemes for time-marching to the steady states is a subject 
of future research. 


29 


Acknowledgment The major portion of the work of the second author was performed as a visiting 
scientist at RIACS, NASA Ames Research Center. Special thanks to Tom Coakley and Marcel 
Vinokur for their critical review of the manuscript. 

References 

Bayliss, A., Parikh, P., Maestrello, L. and Turkel, E. (1985), “A Fourth-Order Scheme for the 
Unsteady Compressible Navier Stokes Equations, ICASE Report No. 85-44. 

Butcher, J.C. (1987), Numerical Analysis of Ordinary Differential Equations , John Wiley & 
Son, Chichester. 

Carpenter, M., Gottlieb, D. and Abarbanel, S. (1993), “The Stability of Numerical Boundary 
Treatments for Compact High-Order Finite-Difference Schemes,” J. Comput. Phys., Vol. 108, pp. 
272-295. 

Carpenter, M. and Kennedy, C.A. (1994), “ Fourth- Order 2N Storage Runge-Kutta Schemes,” 
NASA TM- 109112, June 1 994. 

Carpenter, M. and Casper, J.H. (1997), “The Accuracy of Shock-Capturing in Two Spatial 
Dimensions,” Proceedings of the 13th AIAA CFD Conference, AIAA-97-2107, Snowmass, CO, 
June 29 - July 2, 1997, pp. 488-498. 

Ciment, M. and Leventhal, H. (1975), “Higher Order Compact Implicit Schemes for the Wave 
Equation,” Math. Comp., Vol. 29, No. 132, pp 985-994. 

Coleman, G. N., Kim, J. and Moser, R. D. (1995), “A Numerical Study of Turbulent Supersonic 
Isothermal-Wall Channel Flow,” J. Fluid Mech., Vol. 305, pp. 159-183. 

Davoudzadeh, F., McDonald, H and Thompson, B.E., (1995), “Accuracy Evaluation of 
Unsteady CFD Numerical Schemes by Vortex Preservation,” Computers & Fluids, Vol. 24, No. 
8, pp. 883-895. 

Donat, R. (1994), “Studies on Error Propagation for Certain Nonlinear Approximations to 
Hyperbolic Equations: Discontinuities in Derivatives,” SIAM J. Numer. Anal., Vol. 31, No. 3, 
pp. 655-679. 

Engquist, B. and Sjogreen, B. (1995), “High Order Shock Capturing Methods,” CFD Reviews , 
Hafez & Oshima, Eds., John Wiley, New York, pp. 210-233. 

Fu, D. and Ma, Y. (1997), “A High Order Accurate Difference Scheme for Complex Flow 
Fields,” J. Comput. Phys., vol. 134, pp. 1-15. 

Gear, C.W. (1971), Numerical Initial Value Problems in Ordinary Differential Equations , 
Prentice-Hall, Englewood Cliffs. 

Gerritsen, M. (1996), “Designing an Efficient Solution Strategy for Fluid Flows,” PhD Thesis, 
SCCM Program, Stanford University, Oct. 1996. 

Gottlieb, S. and Shu, C.-W. (1996), “Total Variation Diminishing Runge-Kutta Schemes,” 
ICASE Report No. 96-50. 

Gottlieb D. and Turkel, E., “Dissipative Two-Four Methods for Time Dependent Problems, 
Math. Comp., Vol. 30, pp. 703-723. 



30 


Gustafsson, B. and Olsson, O. (1995), “Fourth-Order Difference Method for Hyperbolic 
IB VPs,” J. Comp. Phys., Vol. 1 177, pp. 300-317. 

Harten, A. (1978), “The Artificial Compression Method for Computation of Shocks and 
Contact Discontinuities: III. Self-Adjusting Hybrid Schemes,’’ Math. Comp., Vol. 32, No. 142, 
pp. 363-389. 

Harten, A. ( 1 98 1 ), ‘ ‘On the Symmetric Form of Systems for Conservation Laws with Entropy,’ ’ 
ICASE Report 81-34, NASA Langley Research Center. 

Harten, A. (1983), “A High Resolution Scheme for Computation of Weak Solutions of 
Hyperbolic Conservation Laws,’’ J. Comp. Phys., Vol. 49, pp. 35-393. 

Harten, A. (1984), “On a Class of High Resolution Total- Variation-Stable Finite-Difference 
Schemes,’’ SIAM J. Num. Anal., Vol. 21, pp. 1-23. 

Harten, A. and Hyman, J.M. (1983), “A Self-Adjusting Grid for the Computation of Weak 
Solutions of Hyperbolic Conservation Laws,” J. Comp. Phys., Vol. 50, pp. 235-269. 

Harten, A. and Osher, S. (1987), “Uniformly High-Order Accurate Nonoscillatory Schemes 
I,” SIAM J. Num. Analy. Vol. 24, No. 2, pp. 279-309. 

Hirsh, R.S. ( 1 975), ‘ ‘Higher Order Accurate Difference Solutions of Fluid Mechanics Problems 
by a Compact Differencing Technique,” J. Comp. Phys., Vol. 19, pp 90-109. 

Hixon, R. and Turkel, E. (1998), “High- Accuracy Compact MacCormack-type Schemes 
for Computational Aeroacoustics,” AIAA paper 98-0365, 36th Aerospace Sciences Meeting & 
Exhibit, Jan. 12-15, Reno, NV. 

Jameson, A., Schmidt W. and Turkel, E. (1981), “Numerical Solutions of the Euler Equations 
by Finite Volume Methods Using Runge-Kutta Time-Stepping Schemes,” AIAA Paper No. 
81-1259. 

Jin, S. and Liu, J.-G. (1996), “The Effects of Numerical Viscosities, I: the Slowly Moving 
Shocks,” J. Comp. Phys., Vol. 126, pp. 373-389. 

Kim, J., Moin, P. and Moser, R.D. ( 1 987), ‘ ‘Turbulence Statistics in Fully Developed Turbulent 
Channel Flow at Low Reynolds Number,” J. Fluid Mech., Vol. 177, pp. 133-166. 

Lambert, J.D. (1973), Computational Methods in Ordinary Differential Equations, John Wiley, 
New York. 

Lee, S., Lele, S. K. and Moin, P. (1997) “Interaction of Isotropic Turbulence with Shock 
Waves: Effect of Shock Strength,” J. Fluid Mech., vol. 340, pp. 225-247. 

Lele, S. (1992), “Compact Finite Difference Schemes with Spectral-Like Resolution,” J. 
Comput. Phys., Vol. 103, pp. 16-42. 

Lumpp, T. (1996a), “A Critical Review of High-Order Accurate Finite Volume ENO- 
Schemes,” First Inti. Symposium on Finite Volumes for Complex Applications, Rouen, France, 
July 14-17. 

Lumpp, T. (1996b), “Compressible Mixing Layer Computations with High-Order ENO 
Schemes,” 15th Inti. Conf. on Num. Meth. in Flui Dynamics, Monterey, June 1996. 

Moin, P. and Mahesh, K. 1998 Direct Numerical Simulation: a Tool in Turbulence Research,” 
Ann. Rev. Fluid Mech., Vol. 30, pp. 539-578. 


31 


Olsson, P. (1995a), “Summation by Parts, Projections and Stability,’’ Math. Comp., Vol. 64, 
pp. 1035-1065. 

Olsson, P. (1995b), “Summation by Parts, Projections and Stability II.,” Math. Comp., Vol. 
64, p212. 

Olsson, P. ( 1 995c), “Summation by Parts, Projections and Stability III.,’ ’ RIACS Report 95-06, 
NASA Ames Research Center. 

Roe, P.L. (1981), “Approximate Riemann Solvers, Parameter Vectors, and Difference 
Schemes,” J. Comp. Phys., Vol. 43, pp. 357-372. 

Sandham, N.D. and Howard, R.J.A. (1997) “Direct Simulation of Turbulence Using Massively 
Parallel Computers,” Parallel Computational Fluid Dynamics ’97, Ed. A. Ecer et al., Elsevier, to 
appear. 

Sandham, N.D. and Yee, H.C. (1989), “A Numerical Study of a Class of TVD Schemes for 
Compressible Mixing Layers,” NASA TM- 102 194, May 1989. 

Sandham, N.D. and Yee, H.C. (1998) “Performance of Low Dissipative High Order TVD 
Schemes for Shock-Turbulence Interactions,” RIACS Technical Report 98.10, April 1998, NASA 
Ames Research Center. 

Shu , S.-W. (1989), “Efficient Implementation of Essentially Non-Oscillatory Shock-Capturing 
Schemes, II,” J. Comput. Phys., Vol. 83, pp. 32-78. 

Shu , S.-W., Erlebacher, G., Zang, T.A., Whitaker, D. and Osher, S. (1991), “High-Order 
ENO Schemes Applied to Two- and Three-Dimensional Compress Flow,” ICASE Report 187562, 
NASA Langley Research Center. 

Shu, S.-W., (1996), Private Communication. 

Strand, B. (1994), “Summation by Parts for Finite Difference Approximations for d/dz,” J. 
Comput. Phys., vol 1 10, pp. 47-67. 

Strang, G. (1968), “On the Construction and Comparison of Difference Schemes,” SIAM J. 
Numer. Anal., Vol. 5, pp. 506-517. 

Toro, E.F. (1991), “Viscous Flux Limiters,” Proceedings of the 9th GAMM-Conference on 
Numerical Methods in Fluid Mechanics,” Notes on Numerical Fluid Mechanics, Vol. 35, Vos, 
J.B, Rizzi, A. and Ryhming, I.L., Editors, pp. 592-600. 

Workshop on Benchmark Problems in Computational Aeroacoustic (1995), NASA Conference 
Publication 3300, edited by Hardin, J.C., Ristorcelli, J.R. and Tam, C.K.W. 

Yee, H.C., Warming, R.F. and Harten, A. (1985), “Implicit Total Variation Diminishing (TVD) 
Schemes for Steady-State Calculations,” J. Comput. Phys., Vol. 57, pp. 327-360. 

Yee, H.C. (1985b), “On Symmetric and Upwind TVD Schemes,” Proceedings of the 6th 
GAMM-Conference on Numerical Methods in Fluid Mechanics, Sept. 1985, Gottingen, Germany. 

Yee, H.C. and Harten, A. (1985), “Implicit TVD Schemes for Hyperbolic Conservation Laws 
in Curvilinear Coordinates,” AIAA Paper No. 85-1513, Proceedings of the AIAA 7th CFD 
Conference, Cinn., Ohio, July, 1985, also AIAA J., Vol. 25, No. 2, 1987, pp. 266-274. 

Yee, H.C. (1986) “Linearized Form of Implicit TVD Schemes for Multidimensional Euler 
and Navier-Stokes Equations,” Computers and Mathematics with Applications, Vol. 12A, pp. 


32 


413-432, 1986. 

Yee, H.C. (1987) “Construction of Explicit and Implicit Symmetric TVD Schemes and Their 
Applications,” J. Comput. Phys., Vol. 68, pp. 151-179, 1987; also NASATM-86775 July 1985. 

Yee, H.C. and Shinn, J. (1989), “Semi-Implicit and Fully Implicit Shock-Capturing Methods 
for Hyperbolic Conservation Laws with Stiff Source Terms,” AIAA-87-1 1 16, AIAA 8th CFD 
Conference, June 9-11, 1987, Hawaii, also, AIAA J., Vol. 27, No. 3, pp. 299-307. 

Yee, H.C. (1989), “A Class of High-Resolution Explicit and Implicit Shock-Capturing 
Methods,” VKI Lecture Series 1989-04 March 6-10, 1989, also NASA TM- 101088, Feb. 1989. 

Yee, H.C., Klopfer, G.H. and Montagne, J.-L. (1990), “High-Resolution Shock-Capturing 
Schemes for Inviscid and Viscous Hypersonic Flows,” J. Comput. Phys., Vol. 88, pp. 31-61. 

Yee, H.C. (1995), “Explicit and Implicit Compact High-Resolution Shock-Capturing Methods 
for Multidimensional Euler Equations I: Formulation,” NASA TM-1 10364, August 1995, also J. 
Comput. Phys., Vol. 131, No. 1, pp. 216-232. 

Young, V.Y.C and Yee, H.C. (1987), “Numerical Simulation of Shock Wave Diffraction by 
TVD Schemes,” AIAA Paper No., 87-01 12. 



Method 

Order 

(Euler) 

Order 

(viscous) 

Shock- 

capturing 

Artificial 

compression 

Compact 

CEN44 

4 

4 

No 

No 

No 

TVD22 

2 

2 

Yes 

No 

No 

TVD44 

4 

4 

Yes 

No 

No 

TVD66 

6 

6 

Yes 

No 

No 

ACM22 

2 

2 

Yes 

Yes 

No 

ACM44 

4 

4 

Yes 

Yes 

No 

ACM66 

6 

6 

Yes 

Yes 

No 

ACM44C 

4 

4 

Yes 

Yes 

Yes 

ACM66C 

6 

6 

Yes 

Yes 

Yes 


Table 3.1. Notation for numerical methods. Order of accuracy refers to the formal order of the 
base scheme. 
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Table 3.2. 


Property 

(i) 

(2) 

( 3 ) 

( 4 ) 

( 5 ) 

u — velocity 

3.0000 

2.0000 

2.9709 

2.9792 

1.9001 

v — velocity 

0.0000 

0.0000 

-0.1367 

-0.1996 

-0.1273 

6 (degrees) 

0.0000 

0.0000 

2.6343 

3.8330 

3.8330 

density p 

1.6374 

0.3626 

2.1101 

1.8823 

0.4173 

pressure p 

0.3327 

0.3327 

0.4754 

0.4051 

0.4051 

sound speed c 

0.5333 

1.1333 

0.5616 

0.5489 

1.1658 

Mach number \M\ 

5.6250 

1.7647 

5.2956 

5.4396 

1.6335 


Flow properties for the shock-wave/shear-layer test case in various regions of the flow: 
(1) upper stream inflow, (2) lower stream inflow, (3) upper stream after oblique shock, 
(4) upper stream after expansion fan, (5) lower stream after shock wave. 
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Property 

[8] 

Small domain 

Centreline velocity 

18.20 

18.23 

Mean velocity 

15.63 

15.69 

Shape factor 

1.62 

1.61 

Skin friction 

0.00818 

0.00811 


Table 3.3. Comparison of incompressible channel flow statistics between large and small 
domain calculations. 
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Fig. 3.1. Density contours 



Density cut at y=0 Density cut at y=0 


Comparison of Various Orders of TVD & ACM Methods 

Stationary Vortex Evolution, 401 x 81 



Fig. 3.2. Stationary Vortex: Comparison of the various orders of TVD and ACM methods with 
the exact solution, illustrated by density profiles at the centerline y = 0, at t = 50 and 
t = 100 for a 401 x 81 grid (k = 0.05). 
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Fig. 3.6. Stationary Vortex: Comparison of the exact solution with ACM66 at t ~ 100 on the 
small computational box (81 x 81) and a larger computational box (401 x 81) (« = 0.05). 



Fig. 3.7. Stationary Vortex: Comparison of the exact solution with ACM44 and ACM66 at 
t = 100 for the small computational box using a 81 x 81 grid (/c = 0.035). 





Fig. 3.8. Stationary Vortex: Comparison of the exact solution with ACM44 using two different 
flux limiters (lim3 - (2.25f); lim5 - (2.25h)) at t = 100 for the small computational box 
using a 81 x 81 grid ( k = 0.05). 







Fig. 3.10. Effect of order of accuracy on TVD methods (TVD22, TVD44 and TVD66), compared 
with the ACM44 solution at t = 160, illustrated by temperature contours at t = 160 for 
a 101 x 101 grid. 
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Fig. 3.12. Effect of k in the ACM44 scheme on the fine scale flow structure, illustrated by 
temperature contours at t = 160: (a) reference solution using n = 0.7 for all fields, (b) 
reference solution using « = 0.7 for the nonlinear fields and k — 0.35 for the linear 
fields, (c) k> = 0.7 for the nonlinear fields and k = 0.35 for the linear fields for a 
101 x 101 gird and (d) /c = 0.7 for the nonlinear fields and k — 0.0 for the linear fields 
for a 41 x 41 grid. 
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Fig. 3.14. Effect of k (with k = 0.0 for the linear fields) in the ACM44 scheme on the fine scale 
flow structure, illustrated by temperature contours at t = 160: (a) n = 0.7 for the 
nonlinear fields, (b) k = 1.0 for the nonlinear fields, and (c) k = 2.0 for the nonlinear 
fields for a 101 x 101 grid. 
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Fig. 3. 15. Effect of flux limiters in the ACM44 scheme on the solution resolution, illustrated by 
temperature contours at t = 160 : (a) liml (2.24d), (b) lim2 (2.24e), (c) lim4 (2.24g), (d) 
lim5 (2.24h) for a 101 x 101 grid with k = 0.7 for the nonlinear fields and k — 0.35 for 
the linear fields. 
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Fig. 3.16. Effect of 6 ( (smu) 2 ) in the ACM44 scheme on the solution resolution, illustrated by 
temperature contours at t - 160: (a) 6 = 1/16, (b) 6 = 0.0225, (c) 6 - 0.01, (d) 6 = 0.0 
for a 101 x 101 grid with k = 0.7 for all fields. 
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Fig. 3.18. Effect of order of accuracy on ACM22/*/ 2 , ACM44A/ 2 and ACM66/*/ 2 on the solution 
resolution, illustrated by temperature contours at t — 160: (a) ACM44 (b) ACM22/i/ 2 
(c) ACM44/*/ 2 (d) ACM66/*/ 2 fora 101 x 101 grid with k = 0.7 for the nonlinear fields 
and k = 0.35 for the linear fields. 
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Fig. 3.21. Comparison of density contours at t - 120 for the shock-shear-layer test case: (a) 
ACM22, (b) ACM44, and (c) ACM66 for a 321 x 81 grid with a = 0.7 for the nonlinear 
fields and k — 0.35 for the linear fields. 








Fig. 3.22. Effect of S ( (smu) J ) in the ACM44 scheme on the solution resolution, illustrated by 
density contours at t = 120 (a) 6 = 0.0225, (b) 6 = 0.01, (c) 6 - 0.0 for a 321 x 81 
grid with k — 0.7 for all fields. 
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Fig. 3.23. Comparison of density contours at t = 120: (a) ACM44/i/ 2 (641 x 161) (b) TVD44/i/ 2 
(321 x 81) (c) ACM44 lu 2 (321 x 81) with « = 0.7 for all fields. 








Fig. 3.24. Compressible turbulent channel flow: (a) mean velocity profile (solid line) compared to 
the linear law u + = y + and the semi-logarithmic law of the wall, with Karman constant 
0.41 and additive constant 5.5, and (b) Favre-averaged (solid line) and total stresses 
(dashed line) compared to the linear result (chain dotted line). 
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